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Abstract. We present a systematic methodology to develop high order accurate numerical 
approaches for linear advection problems. These methods are based on evolving parts of the jet 
of the solution in time, and are thus called jet schemes. Through the tracking of characteristics 
and the use of suitable Hermite interpolations, high order is achieved in an optimally local 
fashion, i.e. the update for the data at any grid point uses information from a single grid 
cell only. We show that jet schemes can be interpreted as ad vect-and— project processes in 
^ function spaces, where the projection step minimizes a stability functional. Furthermore, this 

Q function space framework makes it possible to systematically inherit update rules for the higher 

derivatives from the ODE solver for the characteristics. Jet schemes of orders up to five are 
applied in numerical benchmark tests, and systematically compared with classical WENO finite 
00 difference schemes. It is observed that jet schemes tend to possess a higher accuracy than 

T— H WENO schemes of the same order. 

1. Introduction 

In this paper we consider a class of approaches for Unear advection problems that evolve parts 
of the jet of the solution in time. Therefore, we call them jet sc/iemesj^ As we will show, the idea 
I— I of tracking derivatives in addition to function values yields a systematic approach to devise high 

order accurate numerical schemes that are very localized in space. The results presented are a 
^ generalization of gradient-augmented schemes (introduced in 21 ) to arbitrary order. 

In this paper we specifically consider the linear advection equation 

m 0t + i/-V(/) = o (1) 

for 0, with initial conditions 0) = Furthermore, we assume that the (given) velocity 

field v{x, t) is smooth. Equation ([T]) alone is rarely the central interest of a computational project. 
I However, it frequently occurs as a part of a larger problem. One such example is the movement of 

a front under a given velocity field (in which case (f> would be a level set function '22'). Another 
example is that of advection-reaction, or advection-diffusion, problems solved by fractional steps. 
Here, we focus solely on the advection problem itself, without devoting much attention to the 
background in which it may arise. However, we point out that generally one cannot, at any time, 
simply track back the solution to the initial data. Instead, the problem background generally 
enforces the necessity to advance the approximate solution in small time increments At. 

The accurate (i.e. high-order) approximation of ([T]) on a fixed grid is a non-trivial task. Com- 
monly used approaches can be put in two classes. One class comprises finite difference and finite 
volume schemes, which store function values or cell averages. These methods achieve high order 
accuracy by considering neighborhood information, with wide stencils in each coordinate direc- 
tion. Examples are ENO [25] or WENO [TS] schemes with strong stability preserving (SSP) time 
stepping [lTl[24l[25] or flux limiter approaches |28, 29 . Due to their wide stencils, achieving high 
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^The term "jet scheme" exists in the fields of algebra and algebraic geometry, introduced by Nash 20 and popu- 
larized by Kontsevich |16| . as a concept to understand singularities. Since we are here dealing with a computational 
scheme for a partial differential equation, we expect no possibility for confusion. 
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order accurate approximations near boundaries can be challenging with these methods. In addi- 
tion, the generalization of ENO/WENO methods to adaptive grids (quadtrees, octrees) [H HH] is 
non-trivial. The other class of approaches is comprised by semi-discrete methods, such as discon- 
tinuous Galerkin (DG) [71 [Ml [23]. These achieve high order accuracy by approximating the flux 
through cell boundaries based on high order polynomials in each grid cell. Generally, DG methods 
are based on a weak formulation of ([T]) , and the flow between neighboring cells is determined by a 
numerical flux function. All integrals over cells and cell boundaries are approximated by Gaussian 
quadrature rules, and the time stepping is done by SSP Runge-Kutta schemes [TTl [531 [5S] . While 
the DG formalism can in principle yield any order of accuracy, its implementation requires some 
care in the design of the data structures (choice of polynomial basis, orientation of normal vectors, 
etc.). Furthermore, SSP Runge-Kutta schemes of an order higher than 4 are quite difficult to 
design p[TU]. 

The approaches (jet schemes) considered here fall into the class of semi-Lagrangian approaches: 
they evolve the numerical solution on a fixed Eulerian grid, with an update rule that is based on the 
method of characteristics. Jet schemes share some properties with the methods from both classes 
described above. Like DG methods, they are based on a high order polynomial approximation in 
each grid cell, and store derivative information in addition to point values. However, instead of 
using a weak formulation, numerical flux functions, Gaussian quadrature rules, and SSP schemes, 
here characteristic curves are tracked — which can be done with simple, non-SSP, Runge-Kutta 
methods. Furthermore, unlike semi-discrete methods, jet schemes construct the solution at the 
next time step in a completely local fashion (point by point). Each of the intermediate stages of 
a Runge-Kutta scheme does not require the reconstruction of an approximate intermediate rep- 
resentation for the solution in space. This property is shared with Godunov-type finite volume 
methods. Fundamental differences with finite volume methods are that jet schemes are not conser- 
vative by design, and that higher derivative information is stored, rather than reconstructed from 
cell averages. Finally, like many other semi-Lagrangian approaches, jet schemes treat boundary 
conditions naturally (the distinction between an ingoing and an outgoing characteristic is built 
into the method), and they do not possess a Courant-Friedrichs-Lewy (CFL) condition that re- 
stricts stability. This latter property may be of relevance if the advection equation ([T]) is one step 
in a more complex problem that exhibits a separation of time scales. 

For advection problems ([!]), jet schemes are relatively simple and natural approaches that yield 
high order accuracy with optimal locality: to update the data at a given grid point, information 
from only a single grid cell is used [21] . In addition, the high order polynomial approximation 
admits the representation of certain structures of subgrid size — see [U] for more details. 

Jet schemes are based on the idea of advect-and-project: one time step in the solution of 
Equation ([T]) takes the form 

^PoAt+At.t^, (2) 

where: (i) (p™ denotes the solution at time tm — with tn+i — i„-l-Ai, (ii) At+At, t is an approximate 
advection operator, as obtained by evolving the solution along characteristics using an appropriate 
ODE solver, and (iii) P is a projection operator, based on knowing a specified portion of the jet of 
the solution at each grid point in a fixed cartesian grid. At the end of the time step, the solution 
is represented by an appropriate cell based, polynomial Hermite interpolant — produced by the 
projection P. We use polynomial Hermite interpolants because they have a very useful stabilizing 
property: in each cell, the interpolant is a polynomial minimizer for the norm of a certain high 
derivative of the interpolated function. This controls the growth of the derivatives of the solution 
(e.g.: oscillations), and thus ensures stability. 

Equation ^ is not quite a numerical scheme, since advecting the complete function 0" would 
require a continuum of operations. However, it can be easily made into one. In order to be able to 
apply P, all we need to know is the values of some partial derivatives of At+At. t 0" (the portion 
of the jet that P uses) at the grid points. These can be obtained from the ODE solver as follows: 
consider the formula (provided by using the ODE solver along characteristics) that gives the value 
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of Af^^t tcf)" at any point (in particular, the grid points). Then take the appropriate partial 
derivatives of this formula; this gives an update rule that can be used to obtain the required data 
at the grid points. This process provides an implementable scheme that is fully equivalent to the 
continuum (functional) equation ([2]). We call such a scheme superconsistent, since it maintains 
functional consistency between the function and its partial derivatives. 

Finally, we point out that no special restrictions on the ODE solver used are needed. This follows 
from the minimizing property of the Hermite interpolants mentioned earlier, and superconsistency 
— which guarantees that the property is not lost, as the time update occurs in the functional sense 
of Equation ^ . Thus, for example, regular Runge-Kutta schemes can be used with superconsistent 
jet schemes — unlike WENO schemes, which require SSP ODE solvers to ensure total variation 
diminishing (TVD) stability p] . 

This paper is organized as follows. The polynomial representation of the approximate solution 
is presented in § [2j There we show how, given suitable parts of the jet of a smooth function, a 
cell-based Hermite interpolant can be used to obtain a high order accurate approximation. This 
interpolant gives rise to a projection operator in function spaces, defined by evaluating the jet of 
a function at grid points, and then constructing the piecewise Hermite interpolant. In § |3j the 
jet schemes' advect-and-project approach sketched above is described in detail. In particular, we 
show that superconsistent jet schemes are equivalent to advancing the solution in time using the 
functional equation ^ . This interpretation is then used to systematically inherit update rules for 
the solution's derivatives, from the numerical scheme used for the characteristics. Specific two- 
dimensional schemes of orders 1, 3, and 5 are constructed. These are then investigated numerically 
in § |4] for a benchmark test, and compared with classical WENO schemes of the same orders. 
Boundary conditions and stability are discussed in § |3.4| and § |3.5[ respectively. 



2. Interpolations and Projections 



In this section we discuss the class of projections P that are used and required by jet schemes — 
see Equation We begin, in § 2.1 by presenting the cell-based generalized Hermite interpolation 
problem of arbitrary order in any number of dimensions. In § 2.2 we construct, on an arbitrary 
cartesian grid, global interpolants — using the cell-based generalized Hermite interpolants of § 2.1 
Then we introduce a stability functional, which is the key to the stability of the superconsistent 
jet schemes 

In 



2.3 two different types of portions of the jet of a function (corresponding 
These are the total and the partial A:-jets 



2.4 



Next, in 

to different kinds of projections) are defined 
the global interpolants defined in § 2.2 are used to construct (in appropriate function spaces) the 
projections which are the main purpose of this section. Finally, in § 2.5 we discuss possible ways 
to decrease the number of derivatives needed by the Hermite interpolants, using cell based finite 
differences; the notion of an optimally local projection is introduced there. 



2.1. Cell-Based Generalized Hermite Interpolation. We start with a review of the general- 
ized Hermite interpolation problem in one space dimension. Consider the unit interval [0, 1]. For 
each boundary point q e {0, 1} let a vector of data (0q, 0^, . . . , (j>1) be given, corresponding to all 
the derivatives up to order k of some sufficiently smooth function = 4'{x) — the zeroth order 
derivative is the function itself. Namely 

< = iir^i<l) Vge{0,l}, a e {0,1,..., A:}. 

In the class of polynomials of degree less than or equal to n — 2k + 1, this equation can be used 
to define an interpolation problem with a unique solution. This solution, the n*'' order Hermite 
interpolant, can be written as a linear superposition 

qe{o.,i} ae{Q.,...M 
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of basis functions w'^ ^, each of which solves the interpolation problem 

(s)"'<,a(9')=^a,a"^,,,' Vg'e{0,l}, a' € {0, . . . , fc}, 

where (5 denotes Kroneeker's delta. Hence, each of the 2(fc + 1) = n + 1 basis polynomials equals 1 
on exactly one boundary point and for exactly one derivative, and equals for any other boundary 
point or derivative up to order k. Notice that 

<aW = (-ir<a(l-^)- 

Example 1. The three lowest order cases of generalized Hermite interpolants are given by the 
following basis functions: 

• linear (A; = 0, n = 1): 

w^i.oW = X, 

• cubic (fc = 1, n = 3): 

wI q{x) = - 2x^, 
wl i{x) = —x^ + x^, 

• quintic {k = 2, n = 5): 

^1 = lOx^ - Ibx^ + 6a;^ 
= -Ax^ + 7x'^ -ix^, 

One dimensional Hermite interpolation can be generalized naturally to higher space dimensions 
by using a tensor product approach, as described next. 

In W, consider a p-rectangle (or simply "cell") [ai, 6i] x ••• x [op, bp]. Let Axi = bi — Ui, 
1 < i < p, denote the edge lengths of the p-rectangle, and call h = max^^-^ Axi the resolution. In 
addition, we use the classical multi- index notation. For vectors x & W and a € Nq, define: (i) 
|a| = ELi "i. (ii) = riLi and (iii) = . . . d^" , where 0, = ^. 

Definition 1. A p-n polynomial is a p-variate polynomial of degree less than or equal to n in 
each of the variables. Using multi-index notation, a p-n polynomial can be written as 

with {n+ ly parameters cs- Note that here we will consider only the case where n is odd. 
Example 2. Examples oi p-n polynomials are: 

p = 1 p = 2 p = 3 

• n = l: linear, bi-linear, and tri- linear functions. 

• n = 3: cubic, bi-cubic, and tri-cubic functions. 

• n = 5: quintic, bi-quintic, and tri-quintic functions. 

Now let the p-rectangle's 2^ vertices be indexed by a vector q € {0, 1}^, such that the vertex of 
index q is at position x^- ~ (ai + Axi qi, . . . , + Axp q^). 

Definition 2. For n odd, and a sufficiently smooth function = (/>(a;), the n-data on the p- 
rectangle (defined on the vertices) is the set of (n + 1)^ scalars given by 

0| = a"^(%), (3) 

where qe {0, 1}^ and a e {0, ... , k}P, with k=^. 

Lemma 1. Two p-n polynomials with the same n-data on some p-rectangle, must be equal. 
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Proof. Let (f> be the difference between the two polynomials, which has zero data. We prove that 
= by induction over p. For p = 1, we have a standard 1-D Hermite interpolation problem, 
whose solution is known to be unique. Assume now that the result applies for p — I. In the 
p-rectangle, consider the functions (f>, dp cj), . . . , dp cj) — both at the "bottom" hyperface {xp = Op, 
i.e. Qp = 0) and the "top" hyperface {xp — bp, i.e. qp = 1). For each of these functions, zero 
data is given at all of the corner vertices of the two hyperfaces. Therefore, by the induction 
assumption, all of these functions vanish everywhere on the top and bottom hyperfaces. Consider 
now the "vertical" lines joining a point (xi, . . . , Xp-i) € [oi, 6i] x . . . x [ap_i, in the bottom 

hyperface with its corresponding one on the top hyperface. For each of these lines we can use the 
p = 1 uniqueness result to conclude that = identically on the line. It follows that = 
everywhere. □ 

Theorem 2. For any arbitrary n-data (n odd) on some p-rectangle, there exists exactly one p-n 
polynomial which interpolates the data. 

Proof. Lemma [T] shows that there exists at most one such polynomial. The interpolating p-n 
polynomial is explicitly given by 

E E 4 < 5(^1. (4) 

q<£{0,l}'' (3e{0, fe}J' 

where the 3(2?) are p-n polynomial basis functions that satisfy 
They are given by the tensor products 

i=l 

where = is the relative coordinate in the p-rectangle, and the ^ are the univariate 

basis functions defined earlier. □ 



Next we show that the p-n polynomial interpolant given by Equation Q is a (n + 1)*'* order 
accurate approximation to any sufficiently smooth function (f) it interpolates on a p-rectangle. For 
convenience, we consider a p-cube with Axi = • • • = Axp = h. In this case, the interpolant in Q 
becomes 

E E n<o.te)- (5) 

qe{o,i}p 3e{o, ...,k}P i=l 

Lemma 3. Let the data determining the p-n polynomial Hermite interpolant be known only up to 
some error. Then Equation ([5| yields the interpolation error 

ge{0,l}P Se{Q....,k}P / 

where the notation 5u indicates the error in some quantity u. In particular, if the data (f)'i are 
known with O accuracy, then S (d^Hn) = O 

Theorem 4. Consider a sufficiently smooth function (j), and let 'H,i(a') be the p-n polynomial that 
interpolates the data given by 4> on the vertices of a p-cube of size h. Then, everywhere inside the 
p-cube, one has 

a"-H„-a"0 = o(/i"+i-i"i) , (6) 

where the constant in the error term is controlled by the (n-\- 1)"* derivatives of (p. 
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Proof. Let Q be the degree n polynomial Taylor approximation to </>, centered at some point inside 
the p-cube. Then, by construction: (i) d°'Q — d°'(f> = O In particular, the data for 

Q on the p-cube is related to the data for T-Ln (same as the data for 0) in the manner specified 
in Lemmajs] Thus: (ii) d^Hc — d°''Hn = O where T-Lq is the p-n polynomial that 

interpolates the data given by Q on the vertices of the p-cube. However, from Lemma [ij (iii) 
= since Q is a p-n polynomial. From (i), (ii), and (iii) Equation (|6| follows. □ 

In conclusion, we have shown that: the p-n polynomial Hermite interpolant approximates 
smooth functions with (n -1-1)"* order accuracy, and each level of differentiation lowers the order 
of accuracy by one level. Further: in order to achieve the full order accuracy, the data cfiL must 
be known with accuracy O . 

2.2. Global Interpolant. Consider a rectangular computational domain f2 C in p spacial 
dimensions, equipped with a regular rectangular grid. Assume that at each grid point x^?,, labeled 
by m e ZP, a vector of data is given, where o? € {0, . . . , k}^ — for some k G Nq. Given this 
grid data, we define a global interpolant H : $7 — >■ M, which is a piece-wise p-n polynomial (with 
n = 2fc-|- 1). On each grid cell Ti. is given by the p-n polynomial obtained using Equation Q, with 
q related to the grid index rfi by: q — rh — mo — where toq is the cell vertex with the lowest values 
for each component of to. Note that H is C°° inside each grid cell, and C'^ across cell boundaries. 
However, in general, H is not C'^^^. 

Remark 1. The smoothness of H is biased in the coordinate directions. All the derivatives that 
appear in the data vectors, i.e. d^H for a e {0, ... , k}P, are defined everywhere in il. Further- 
more, at the grid points, d°"H{x^) = (j)^. In particular, all the partial derivatives up to order k 
are defined, and continuous. However, not all the partial derivatives of orders larger than k are 
defined. In general d" H, with a. S {0, . . . , k -\- 1}^, is piece-wise smooth — with simple jump 
discontinuities across the grid hyperplanes which are perpendicular to any coordinate direction xe 
such that ag ^ k + 1. Derivatives 9" H, where > fc -|- 1 for some 1 < ^ < p, generally exist only 
in the sense of distributions. 

Definition 3. Any (sufficiently smooth) function (/) : J7 — ^ K defines a global interpolant as 
follows: at each grid point x^, evaluate the derivatives of cj), as by Definition [2j to produce a data 
vector. Namely: 0^ = (3"0(x^i) VcJ G {0, . . . , k}P. Then, use these values as data to define Hif, 
everywhere. 

Definition 4. For any sufficiently smooth function 4>, define the stability functional by 

m= J^{dUix)y dx, (7) 

where f3 = /3{k, p) is the p-vector /3 ^ {k + 1, . . . ,k + l). Of course, J- also depends on fc e No and 
r2 C M^, but (to simplify the notation) we do not display these dependencies. 

Theorem 5. Replacing a sufficiently smooth function (j) by the interpolant T-L^ does not increase 
the stability functional: ^[H^] < J^[<p] ■ In fact: T-L^ minimizes T , subject to the constraints given 
by the data d" 4>{x^), and the requirement that the minimizer should be smooth in each grid cell. 

Remark 2. The minimizer is not unique if p > 1. To see this, let fj{xj), 1 < j < p he nonzero 
smooth functions, with fj and all its derivatives vanishing at the the grid points. Then ip = 
^0 + fji^j) 7^ ^</' the same data as "H^ and 0, and J^[ip] — J^iT-Lff,] if p > 1. 

Proof. Note that d^H^ exists and it is piece- wise smooth (it may have simple discontinuities 
across cell boundaries — see Remark [ij). Hence TIT-Lrf,] is defined. Clearly, it is enough to show 
that minimizes in each grid cell Q. Define (p = (p — H^. Since and have the same data 
at the vertices of Q, if has zero data at the vertices. Thus 

Ip = (a%4x)) i^dP^ix)'^ dx = 0, (8) 
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as shown in Lemma m From this equahty it follows that = J-lV.^] + J-if]- Now, since 

^['■p] ^ 0: J' [Hip] < For any other 0* satisfying the theorem statement's constraints for the 

minimizing class, — "H^.. Hence = J-lV.^'] < ^[0*]- □ 



Corollary 6. — is an orthogonal projection with respect to the positive semi-definite qua- 
dratic form associated with J- — namely: the one used by Equation ^ . 

The reason for the name "stability functional" , and the relevance of Theorem [5) will become 
clear later in § 3.5 where the inequality J-lT-L^] < ^[0] is shown to play a crucial role in ensuring 
the stability of superconsistent jet schemes. Theorem [5] states that the Hermite interpolant is 
within the class of the least oscillatory functions that matches the data — where "oscillatory" is 
measured by the functional J-. 



Lemma 7. The equality in ^ appl: 



les. 



Proof. The proof is by induction over p. Without loss of generality, assume Q ~ [0, 1]^ is the unit 
p-cube. For p — 1 the result holds, since fc + 1 integrations by parts can be used to obtain: 

Jo 

where there are no boundary contributions because the data for tp vanishes, and we have used 
that di^~^'^'Htfi{x) = 0. Assume now that the result is true for p — 1 > 0, and do A: + 1 integrations 
by parts over the variable Xp. This yields 



Ip = (-1) 




where BTC stands for "Boundary Terms Contributions", f3' = + . . . , fc + 1) is a {p— l)-vector, 
and we have used that dp^^^T-L^{x) = 0. Furthermore the BTC are a sum over terms of the form 

Ip-i,j,g = I (d^'dl+=+^ (df^'dl'^ ^(f)) Ax, 0<J<k, 0<q<l, 

J Qq 

where Qq stands for the {p — l)-cube obtained from Q by setting either Xp = {q ~ ()) oi Xp — 1 
{q — 1). Now: the data for tp is, precisely, the union of the data for {d^~^ v}^=o in both Qo and 
Qi. Furthermore dp~^^^^ H^, restricted to either Qq or Qii is a (p — l)-n polynomial. Hence we 
can use the induction hypothesis to conclude that Ip-i,j,q = 0, for all choices of j and q. Thus 
Ip ~ BTC = 0, which concludes the inductive proof. □ 



2.3. Total and Partial fc-Jets. The jet of a smooth function (p : fl M is the collection of all 
the derivatives d"(j), where a £ Nq. The interpolants defined in § 2.2 are based on parts of the 



jet of a function, evaluated at grid points. Next we introduce two different notions characterizing 
parts of jet. 

Definition 5. The total k-jet is the collection of all the derivatives where |c5| < k. 

Definition 6. The partial k-jet is the collection of all derivatives where a £ {0, ... , k}^. 



Thus the total fc-jet contains all the derivatives up to the order fc, while the partial fc-jet contains 
all the derivatives for which the partial derivatives with respect to each variable are at most order 
fc. For any given fc > 0, the total fc-jet is contained within the partial fc-jet (strictly if fc > 0), 
while the partial fc-jet is a subset of the total n-jet, with n = pk. 
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2.4. Projections in Function Spaces. The aim of this subsection is to use the global interpolant 
introduced in §|2.2|to define projections of functions — which are needed for the projection step in 



the general class of advect-and-project methods specified in § 3.2 This requires the introduction 
of appropriate spaces where the projections operate. Further, because of the advection that occurs 
between projection steps in the § |3.2| methods, it is necessary that these spaces be invariant under 
diffeomorphisms. Unfortunately, the coordinate bias (see Remark [l]) that the generalized Hermits 
interpolants exhibit causes difficulties with this last requirement, as explained next. 

The functions that we are interested in projecting are smooth advections, over one time step, of 
some global interpolant %. Let V' be such a function. From Remark [l] it follows that is C*^ and 
piece-wise smooth, with singularity hypersurfaces given by the advection of the grid hyperplanes. 
Inside each of the regions that result from advecting a single grid cell, ip is smooth all the way to 
the boundary of the region. Furthermore, at a singularity hypersurface: 

(i) Partial derivatives of -0 involving differentiation normal to the hypersurface of order k or 
lower are defined and continuous. 

(ii) Partial derivatives of ijj involving differentiation normal to the hypersurface of order A: + 1 
are defined, but (generally) have a simple jump discontinuity. 

(iii) Partial derivatives of "0 involving differentiation normal to the hypersurface of order higher 
than k + 1 are (generally) not defined as functions. However, they are defined on each 
side, and are continuous up to the hypersurface. 

Imagine now a situation where a grid point is on one of the singularity hypersurfaces. Then the 
interpolant H^, as given by Definition [3j may not exist — because one, or more, of the partial 
derivatives of ■0 needed at the grid point does not exist. This is a serious problem for the advect- 
and-project strategy advocated in § |3.2[ This leads us to the considerations below, and a recasting 
of Definition [3] (i.e.: Definition [?]) that resolves the problem. 

We are interested in solutions to Equation ([I]) that are are sufficiently smooth. Let then ■0 be 
an approximation to such a solution. In this case, from Theorem |4] we should add to (i - iii) 
above the following 

(iv) The jumps in any partial derivative of (across a singularity hypersurface) cannot be 
larger than 0(/i"+^^''), where n = 2fc + 1 and s is the order of the derivative. 

Then, from Lemma[3j we sec that these jumps are below the numerical resolution — hence, essen- 
tially, not present from the numerical point of view. This motivates the following generalization 
of Definition [3l 

Definition 7. For functions V' that satisfy the restrictions in items (i - iv) above, define the 
global interpolant Ti,^ as in Definition [3] with one exception: whenever a partial derivative 0^ = 
d°'tjj{xrTi) which is needed for the data at a grid point does not exist (because the grid point is on 
a singularity hypersurface), supply a value using the formula 

d"Tp{x^) = i ( limsup5°0(f) + liniinf 9"0(x) ) . (9) 

Remark 3. The value in equation ([9| for 9"0(a:,n) ignores the "distribution component" of 
d°'ipi^m) — i-G. Dirac's delta functions and derivatives, with support on the singularity hypersur- 
face. These components originate purely because of approximation errors: small mismatches in the 
derivatives of the Hermite interpolants across the grid hyperplanes. Thus we expect Equation ^ 
to provide an accurate approximation to the corresponding partial derivative of the smooth func- 
tion that ijj approximates. Note that, while we are unable to supply a rigorous justification for 
this argument, the numerical experiments that we have conducted provide a validation — see §|4j 

Remark 4. An alternative to Definition [7j that also allows the construction of global interpolants 



from advected global interpolants, is introduced in Definition 12 below. 
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Motivated by the discussion above, we now consider two types of projections (and corresponding 
function spaces). Note that in the theoretical formulation that follows below, the "small jumps 
in the partial derivatives" property is not built into the definitions of the spaces. This is not 
untypical for numerical methods. For instance, finite difference discretizations are only meaningful 
for sufficiently smooth functions. However, the finite differences themselves are defined for any 
function that lives on the grid, even though the notion of smoothness does not make sense for such 
grid-functions. 

Definition 8. Let S"" denote the space of of all the L°° functions : ^ M which are v-t\mes 
differentiable a.e., with derivatives in L°° . 

Here differentiable is meant in the classical sense: ^(x + u) = ^j{x) + [D tp) [u\ + ^ (_D^ -0) [u, u] + 
. . . + ^ {D" -0) [u, . . . , u] +o (||m||'^) near a point where the function is differentiable, where (£>'' ip), 
1 < /i < J/, is a /i-linear form. 

Definition 9. Let S''' + be defined by S''' + = S" n , with v = pk. 

Definition 10. For any function £ S'^, and any d? G Nq with \d\ < ly, define d"ip at every point 
Xffi_ in the rectangular grid, via 

d°'i/j{xfi%) — I- I esslimsup9"V'(2?) + essliminf (9"'0(a;) I . (10) 

where ess lim sup and ess lim inf denote the essential upper and lower limits [2] , respectively. Notice 
that this last formula differs from (|9| by the use of essential limits only. 

The definition in Equation ( |10[ ) is not the only available option. One could use the upper limit 
only, or the lower limit only, or some intermediate value between these two — since numerically the 
various alternatives lead to differences that are of the order of the truncation errors. In particular, 
a single sided evaluation is more efficient, so this is what we use in our numerical implementations. 



In what follows we assume that a choice has been made, e.g.: the one given by Equation (10), or 
the one given in Remark [9] 

Definition 11. In S''^'"'", define the projection Pj^ : 5''^'+ — > S^''^ as the application of the 



interpolant defined in § 2.2 i.e. P^ip — "H^ 



Definition 12. In S , define the projection Pk : S — > S by the following steps. 

(1) Given a function ip g S'^, evaluate the total /c-jet at the grid points x,ri- 

(2) Approximate the remaining derivatives in the partial fc-jet by a process like the one described 
in § |2.5[ Notice that, if the approximation of the derivatives of order ^ > fc is done with 
0(ft,"+-'^~^) errors, then (by Lemma^ the full accuracy of the projection is preserved. 



(3) Define Pkip as the global interpolant (see § 2.2) based on this approximate partial fc-jet 



2.5. Construction of the Partial fc-jet from the Total fc-jet. For the advect-and-project 
methods introduced in § |3.2[ the advection of the solution's derivatives is a significant part of the 
computational cost. Furthermore, the higher the derivative, the higher the cost. For example. 



when using P^ with k = 2 and the approach in § 3.3.2 the cost of obtaining the \a\ > k 
derivatives in the partial fc-jet is (at the lowest) about three times that of obtaining the total 
fc-jet. Further: the cost ratio grows, roughly, linearly with fc. Thus it is tempting to design more 
efficient approaches that are based on advecting the total fc-jet only, and then recovering the higher 
derivatives in the partial fc-jet by a finite difference approximation of the total fc-jet data at grid 
points. We call such approaches grid-based finite differences, in contrast to the e- finite differences 
introduced in § |3.3.2[ 

Unfortunately, jet schemes based on grid-based finite difference reconstructions of the higher 
derivatives in the partial fc-jet have one major drawback: the nice stabilizing properties that 



the Hermite interpolants possess (see Theorem [s] and § 3.5) are lost. Hence the stability of 



these schemes is not ensured. Nevertheless, we explore this idea numerically (see § 4.2 for some 
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examples). As expected, in the absence of a theoretical underpinning that would allow us to 
distinguish stable schemes from their unstable counterparts, many approaches that are based on 
grid-based finite differences turn out to yield unstable schemes. However 

• In the case fc = 1, a stable scheme using a grid-based finite difference reconstruction of the 
partial fc-jet, can be given. It is described in Example [3j The particular reconstruction 
used has some interesting properties, which motivate Definition |13| below. 

• In unstable versions of schemes that are based on grid-based reconstructions of higher 
derivatives, the instabilities are, in general, observed to be very weak: even for fairly 
small grid sizes, the growth rate of the instabilities is rather small — an example of this 



phenomenon is shown in § 4.2 Thus it may be possible to stabilize these schemes (e.g. by 



some form of jet scheme artificial viscosity) without seriously diminishing their accuracy. 

Definition 13. We say that a projection is optimally local if, at each grid node, the recovered 
data (e.g. derivatives with |a| > A; in the partial fc-jet) depends solely on the known data (e.g. the 
total fc-jet) at the vertices of a single grid cell. 



Of course, the projections in Definition 11 are by default optimally local (since there is no data 
to be recovered). However, as pointed out above, the idea of optimally local projections that are 
based on the total fc-jet is that they produce more efficient jet schemes than using . In Examplejs] 
below we present an optimally local projection, which corresponds to Pi in Definition |12[ 

Remark 5. An important question is: why is optimal locality desirable? The reason is that 
optimally local projections do not involve communication between cells. This has, at least, three 
advantageous consequences. First, near boundaries, a local formulation simplifies the enforcement 



of boundary conditions (e.g. see § 3.4). Second, in situations where adaptive grids are used, a 
local formulation can make the implementation considerably simpler than it would be for non- 
local alternatives. And third, locality is a desirable property for parallel implementations since 
communication boundaries between processors are reduced to lines of single cells. 

Example 3. In two space dimensions, we can define an optimally local version of the projection 
Pi, as follows below. In this projection we presume that the total 1-jet is given at each grid point. 
To simplify the description, we work in the cell Q = [0, h] x [0, h], with q S {0, 1}^ corresponding 
to the vertex x^= qh. Thus tp'^ indicates the value ipix^). We also use the convention that when 
qi — ^, the quantity is defined or evaluated at Xi — ^ , for i G {1, 2}. 

Define Pi as follows (this is the projection introduced, and implemented, in [21'). To construct 
the bi-cubic interpolant, at each vertex of Q the derivative must be approximated with second 
order accuracy, using the given data at the vertices. To do so: (i) Use centered differences to write 
second order approximations to at the midpoints of each of the edges of Q. For example: 

"ipxy^^ = ji i^x''^^ — V'i^'*'^) and V'iy'"^ = ji — (ii) Use the four values obtained in 

(i) to define a bi-linear (relative to the rotated coordinate system x + y and x — y) approximation 
to ipxy (iii) Evaluate the bi-linear approximation obtained in (ii) at the vertices of Q. Three 
dimensional versions of these formulas are straightforward, albeit somewhat cumbersome. 

We do not know whether stable analogs of Pi exist for P^ with fc > 1. It is plausible that 
a construction similar to the one given in Example |3] (i.e. based on using lower order Hermite 
interpolants in rotated coordinate systems) may yield the desired stability properties. This is the 
subject of current research. 

Remark 6. The optimally local projection Pi with "single cell approximations" uses, at each 
grid point, values for '>p.j;y that are different for each of the adjoining cells. This means that 
(generally) Piip ^ — though Piip G always. However, for sufficiently smooth functions -0, 
the discontinuities in the derivatives of Pi ip across the grid lines cannot involve jumps larger than 
0(/i '*"''), where s is the order of the derivative. This is a small variation of the situation discussed 
from the beginning of § 2.4 through Remark [Sj The same arguments made there apply here. 
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Remark 7. The lack of smoothness of the optimally local projection Pi (i.e. Remark |6]) could be 
eliminated as follows. At each grid node consider all the possible values obtained for ip^y by cell 
based approximations (one per cell adjoining the node). Then assign to ipxy the average of these 
values. This eliminates the multiple values, and yields an interpolant Piip that belongs to C^. 
Formally this creates a dependence of the partial 1-jet on the total 1-jet that is not optimally local. 
However, all the operations done are solely cell-based (i.e.: the reconstruction of the partial 1-jet) 
or vertex based (i.e.: the averaging). Hence, such an approach can easily be applied at boundaries 
or in adaptive meshes. 



3. Advection and Update in Time 



The characteristic form of Equation ([T]) is 

dx 



«(x, t), (11) 



Let X{x, T, t) denote the solution of the ordinary differential equation (111 at time t, when starting 
with initial conditions x at time t = t. Hence X is defined by 

—X{x, T, t) = v{X{x, T, t), t), with X{x, T, t) — X. 

Due to (12), the solution of the partial differential equation ([T]) satisfies 

(/)(f, t + At)^ (t){X{x, t + At, t), t). 

That is: the solution at time t + At and position x is found by tracking the corresponding 
characteristic curve, given by (111, backwards to time t, and evaluating the solution at time t 
there. Introduce now the solution operator S'^-,*, which maps the solution at time t to the solution 
at time r. Then Sr,t acts on a function g{x) as follows 

[Sr,t9){x)^g{X{x,T,t)). (13) 

Hence, the solution of (jlj satisfies 

0(x, t + At) = St+At,t<l>{^^ 

In particular, the solution at time t = n At can be obtained by applying successive advection steps 
to the initial conditions 

(t){x, t) = St^ t-^At o • • • o S2At, At o SAt, *(^)- 

3.1. Approximate Advection. Let X be an approximation to the exact solution X, as arising 
from a numerical ODE solver — e.g. a high order Runge-Kutta method. By analogy to the true 
advection operator (13), introduce an approximate advection operator, defined as acting on a 
function g{x) as follows 

iAr,tg){x)^g{X{x,T,t)). (14) 
The approximate characteristic curves, given by X, motivate the following definition. 

Definition 14. The point Xfoot = X{x, t + At, t) is called the foot of the (approximate) charac- 
teristic through X at time t + At. Note that Xfoot — 2^foot(a?: t, At), but we do not display these 
dependencies when obvious. 



Let now X{x, t + At, t) represent a single ODE solver step for ( 11 1, from t + At to t, starting 
from X — i.e. let At be the ODE solver time step. Then the successive application of approximate 
advection operators 

At, t-At o • • • o A2At, At o AAt, ^(x) (15) 
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yields an approximation to the solution of ([T]) at time t = nAt. In principle this formula provides 
a way to (approximately) solve (11) by taking (for each point x at which the solution is desired 
at time t — n At) n ODE solver steps from time t back to time 0, and then evaluating the initial 
conditions at the resulting position. However, as described in §[l] we are interested in approaches 
that allow access to the solution at each time step (represented on a numerical grid by a finite 
amount of data), and then advance it forward to the next time step. Hence the method provided 



by expression ( 15 ) is not adequate 



3.2. Advect— and— Project Approach. In this approach, an appropriate projection is applied 
at the end of every time step, after the advection. The projection allows the representation of 
the solution, at each of a discrete set of times, with a finite amount of data. Here, we consider 



projections based on function and derivative evaluations at the grid points, as in § 2.4 Thus, let P 
denote a projection operator — e.g. see Definitions [TTj [T2| or[T3| Then the approximate solution 
method is defined by 

0approx(:?,O = {PoAt.t-At)o---o{PoA2At.At)o{PoAAt,o) H^) . (16) 

Namely, to advance from time t to t + At, apply i^o^t+At.t to the (approximate) solution at time 
t. The key simplification introduced by adding the projection step is that: in order to define the 
(approximate) solution at time t + At, only the data at the grid points is needed. At this point, 
all that is missing to make this approach into a computational scheme, is a method to obtain the 



grid point data at time t + At from the (approximate) solution at time t. This is done in § 3.3 



As shown in § 2.1 the P^. projection is an 0(/i""'"^) accurate approximation for sufficiently 
smooth functions, where n = 2fc + 1 and h = max^ Axi. Thus, with At cx h the use of a locally 
(n + 1)^* order time stepping scheme ensures that the full accuracy is preserved. As usual, the 
global error is in general one order less accurate, since 0(1/ h) time steps are required to reach a 
given final time. Hence, a fc-jet scheme can be up to n*'' order accurate. 



3.3. Evolution of the fc-Jet. The projections defined in § 2A^ require knowledge of parts of the 



jet of the (approximate) solution at the grid points. Specifically, P^ (see Definition |ll[) requires 



the partial /c-jet, and Pk (see Definition 12) requires the total fc-jet. 



A natural way to find the jet at time t + At is to consider the approximate advection operator 



(14 1. Since it defines an approximate solution everywhere, it also defines the solution's spacial 



derivatives. Thus, by differentiating (14), update rules for all the elements of the fc-jet can be 



systematically inherited from the numerical scheme for the characteristics. 

Definition 15. Jet schemes for which the update rule for the whole (total or partial) fc-jet is 
derived from one single approximate advection scheme are called superconsistent. 

Remark 8. Non-superconsistent jet schemes can be constructed, by first writing an evolution equa- 
tion along the characteristics for each of the relevant derivatives — by differentiating Equation ([T]), 
and then applying some approximation scheme to each equation. On the one hand, these schemes 
can be less costly than superconsistent ones, because an t*'^ order derivative only needs an accu- 
racy that is £ orders below that of the function value — see Lemma [3) On the other hand, the 
benefits of an interpretation in function spaces are lost, such as the optimal coherence between all 
the entries in the fc-jet, and the stability arguments of § |3.5[ 



Notice that a superconsistent scheme has a very special property: even though it is a fully 
discrete process, it updates the solution in time in a way that is equivalent to carrying out the 



process given by Equation ( 16 ) in the function space where the projection is defined. Below we 
present two possible ways to find the update rule for the fc-jet: analytical differentiation and e- 
finite differences. Up to a small approximation error, both approaches are equivalent. However, 
they can make a difference with the ease of implementation. 
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3.3.1. Analytical differentiation. Let T-l denote the approximate solution at time t. Since each 
time step ends with a projection step, T-L is a piecewise p-n polynomial, defined by the data at 
time t. Using the short notation X = X{x, t + At, t), update formulas for the derivatives (here 
up to second order) are given by: 

At+AtA<l>{S, t) = H{X, t), 

d dX 

T^A+At.t-^lx, t) = — -DHiX, t), 

, . d^X . ldX\ . dX 



-A,+AM -^(5^, t) = inn^ -Dnix, t)+\ — \ ■ D'H{X, t) 



oXiOXj oxi oxj \ oxi I oxj 

where DH is the Jacobian and D^H the Hessian of It should be clear how to continue this 
pattern for higher derivatives. Since "H is a p-n polynomial, its derivatives are easy to compute 
analytically. The partial derivatives of X follow from the ODE solver formulas, as the following 
example (using a Runge-Kutta scheme) illustrates. 

Example 4. When tracking the total 1-jet (this would be called a gradient- augmented scheme, see 
[21]), the Hermite interpolant is fourth order accurate. Hence, in order to achieve full accuracy 
when scaling At cx h, the advection operator should be approximated with a locally fourth order 
accurate time-stepping scheme. An example is the Shu-Osher scheme j25| . which in superconsistent 
form looks as follows. 

xi = X — At V [x, t + At), 
Vfi = I - AtV v{x, t + At), 



3 



ifi - \Atv{xi, t), 



X2 ^ X I 4 

Vf2 = |/ + i Vfi - i At Vfi • Vu(fi, t), 
x+jx2- I At w (fa, i+ 5 At), 



-^foot 



V ffoot = I / + I V X2 - I At V f2 • V t/(f2, i + i At), 
(t){x, t + At) = ■H(xfoot, t), 
(V(/))(f, t + At) = Vffoot • V-H(ffoot, t). 



Remark 9. When Xfoot is on a grid hyperplane, Equation ( 10 ) in Definition 10 would require several 
ODE solves for derivatives that are discontinuous across the hyperplane (one for each of the local 
Hermite interpolants used in the adjacent grid cells). This would augment the computational 
cost without any increase in accuracy. Hence, in our implementations we do only one update, as 
follows: (i) A convention is selected that assigns the points on the grid hyperplanes as belonging 
to (a single) one of the cells that they are adjacent to. (ii) This convention then uniquely defines 
which Hermite interpolant should be used whenever Xfoot is at a grid hyperplane. 

Finally, we point out that 

• The SSP property p2] is not required for the tracking of the characteristics, since this is 
not a process in which TVD stability plays any role. The use of the Shu-Osher scheme [53] 
in Example [4] is solely to illustrate the analytical differentiation process; any Runge-Kutta 
scheme of the appropriate order will do the job. In particular, when using a 5*'' order 
Runge-Kutta scheme (e.g. the 5*'' order component of the Cash-Karp method [5 ) yields 
a globally 5*'' order jet scheme. 

• Jet schemes that result from using Pq^ (p-linear interpolants) do not need any updates of 
derivatives, since they are based on function values only. In the one dimensional case, the 
CIR method by Courant, Isaacson, and Rees is recovered. For these schemes a simple 
forward Euler ODE solver is sufficient. 

• Using analytical differentiation, the update of a single partial derivative is almost as costly 
as updating all the partial derivatives of the same order. In particular, in order to update 
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the partial fc-jet all the derivatives of order less than or equal to pk have to be updated — 
even though only a few with order greater than k are needed. Therefore, it is worthwhile 
to employ alternative ways to approximate the |q?| > k derivatives in the partial fc-jet. 



One strategy is to use the grid-based finite differences described in § 2.5 assuming that a 
stable version can be found. Another strategy, which is essentially equivalent to analytical 
differentiation, is presented next. 

3.3.2. e-finite differences. Here we present a way to update derivatives in a superconsistent fash- 
ion that avoids the problem discussed in the third item at the end of § |3.3.1[ In e-finite differences 
the definition of a superconsistent scheme is applied by directly differentiating the approximately 
advected function At+At. t 0(2?, t), or its derivatives. This is done using finite difference formu- 
las, with a separation e h chosen so that maximum accuracy is obtained. As shown below, 
£-finite differences can be either employed singly, or in combination with analytic differentiation. 
In particular, in order to update the partial fc-jet, one could update the total fc-jet using analyt- 
ical differentiation, and use e-finite differences to obtain the remaining partial derivatives. The 
approach is best illustrated with a few examples. 

Example 5. Here we present examples of how to implement the advect-and-project approach 



given by Equation (16) in two dimensions, with P = or P = P2 (sec Definition 11), using a 
combination of analytical differentiation and e-finite differences. The error analysis below estimates 
the difference between the derivatives obtained using e-finite differences and the values required 
by superconsistency. This difference can be made much smaller than the numerical scheme's 



approximation errors, so that stability (see § 3.5 ) is preserved. In this analysis, S > characterizes 



1,1) + 






M) + 


0(1-1) - 


0(-l-l)) , 


14) _ 


0(1'-1) - 




1,1) _ 


rA(l'-l) 4- 


,A(-i-i)^ . 



the accuracy of the floating point operations. Namely, S is the smallest positive number such that 
1 + (5 7^ 1 in floating point arithmetic. 

• Advect-and-project using P^ . At time t + At and at every grid point {x^, Vm), the 
values for (0, D0, (p^y) must be computed from the solution at time t. To do so, use the 
approximate advection solver X to compute the approximately advected solution (j)'' at the 
four points {x^ + Qi £, Vm + 12 £)^ where q G {—1, 1}^. From these values ^f^^i', 

and 0(~i'~i) we obtain the desired derivatives at the grid point using the O(e^) 
accurate stencils 

0. = ^ (0(^'^^ - 

• Advect-and-project using P2 ■ At time t + At and at every grid point {x, y) = {x„i, j/m), 
the values for (0, D(f>, D^(f>, 4>xxyi 4>xyyi (f'xxyy) must be computed from the solution at 
time t. To do so, use the approximate advection solver X with analytic differentiation, to 
compute (0, D(j), D'^cj)) at the three points (x,„, Um) and {Xm ± Um)- From this obtain 
{4'xxy, 4>xyy, 4>xxyy) at the grid point using O(e^) accurate centered differences in x. 

Error analysis. In both cases P^ and P^ , the errors in the approximations are: £ = O(e^) + 0{S) 
for averages, £ — O(e^) + 0[5/e) for flrst order centered c-differences, and £ = O(e^) -I- 0{5/e^) 
for second order centered e-differences. Thus the optimal choice for e is e = 0{5^/'^)^ which yields 
£ = 0((5i/^) for all approximations. 

Remark 10. In Remark [9] the scenario of a characteristic's foot point falling on a cell boundary is 
addressed. For e-flnite differences the analogous situation is more critical, since several character- 
istics, separated by 0(e) distances, are used. Whenever the foot points for these characteristics 
fall in different cells, the interpolant corresponding to one and the same cell must be used in the 
e-flnite differences. This is always possible, since the interpolants are p-n-polynomials, and thus 
defined outside the cells. 
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Finally, we point out that, when computing derivatives with e-finite differences, round-off error 
becomes important — the more so the higher the derivative. The error for an s*'* order derivative, 
with an order 0{e^) accurate e-finite difference, has the form £ = + 0{6/e''). The best 

choice is then e — 0(6^^'^^'^'^^), which yields £ = 0((5''/*^*+''''). Hence, under some circumstances 
it may be advisable to use high precision arithmetic and thus make S very small. Given their 
relatively simple implementation, e-finite differences can be the best choice in many situations. 

3.4. Boundary Conditions. In this section we illustrate the fact that, due to their use of the 
method of characteristics, jet schemes treat boundary conditions naturally. For simplicity, we 
consider the two dimensional computational domain fl — [0, 1]^, equipped with a regular square 
grid of resolution h = — i.e. the grid points are = hrh, where m G {0, . . . , N}^. In we 
consider the advection equation 

4>t + + V 4>y = 0, with initial conditions 4>{x, y, 0) = '^{x, y). (17) 

Let us assume that u > on both the left and right boundaries, ?; < on the bottom boundary, 
and u > on the top boundary. Hence the characteristics enter the square domain of computation 
only through the left boundary, where boundary conditions are needed. For example: Dirichlct 
(/)(0, y,t) = ^(y, t) or Neumann 0£c(O, y,t) — (^{y, t). We will also consider the special situation 
where u = on the left boundary, in which case no boundary conditions are needed. 



We now approximate (17 1 with a jet scheme, with the time step restricted by At < /i/A, where 



A is an upper bound for + 1'^ — so that for each grid point ||xfoot ~ 2^11 < h. Then Xfoot S ^ 
for any node {Xm, yn) with m ^ 0. Hence: at these nodes the data defining the solution can 
be updated by the process described earlier in this paper. We only need to worry about how to 
determine the data on the nodes for which a; = 0. There are three cases to describe. 

(i) Dirichlct boundary conditions. Then, on x = we have: (j)t — £.t and cjjy — S,y, while (f>x 
follows from the equation: (j)^ = —{(j3t + v(j)y)/u. Formulas for the higher order derivatives 
of (j) at the left boundary can be obtained by differentiating the equation. For example, 
from 

Ipty +Uy(f)x+U (f>xy +Vy(j)y+V <j)yy = 

we can obtain ipxy 

(ii) Neumann boundary conditions. Then, on x — Q, (j) satisfies (j)t -\- v(f)y — — u^. This is 
a lower dimensional advection problem, which can be solvecj^ to obtain ^ = 0(0, y,t). 
Higher order derivatives follow in a similar fashion. 

(iii) M = on the left boundary. Then Xfoot G ^ for the nodes with a; = 0, so that the data 
defining the solution can be updated in the same way as for the nodes with a: > 0. 

3.5. Advection and Function Spaces: Stability. In this section we show how the interpreta- 
tion of jet schemes as a process of advect-and-project in function spaces, together with the fact 
that Hermite interpolants minimize the stability functional in Equation ([t]), lead to stability for jet 
schemes. We provide a rigorous proof of stability in one space dimension for constant coefficients. 
Furthermore, we outline (i) how the arguments could be extended to higher space dimensions, and 
(ii) how the difficulties in the variable coefficients case could be overcome. Note that, even though 
no rigorous stability proofs in higher space dimensions are given here, the numerical tests shown 
in §|4] indicate that the presented jet schemes are stable in the 2-D variable coefficients case. 

3.5.1. Control over averages in 1-D Hermite interpolants. Let k € No, with n = 2k -\- 1. Let T-L be 
an arbitrary l-n-polynomial, and define 



rk^rk{x)^-^—x''+^{l-xf+^ and - ^fe(x) = rf (18) 



"^Noto that jet schemes can be easily generahzed to problems with a source term. That is, replace Equation ([TJ 
<j>t + V ■ V0 = S, where S is a known function. 
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where /^^' indicates the j*^ derivative of a function / = fix). Integration by parts then yields 



{x)H^''+^\x)dx^{^l)^+' / rk{x)n^''+^Hx)dx^O, 



(19) 



where the boundary contributions vanish because of the {k + 1)** order zeros of at a; = and 
at X = 1. Integration by parts also shows that, for any sufficiently smooth function ip = ipix) 



x=l 



ijj{x) dx, 



where we have used that ^j!^^^^^ = (—1)'^+^. This can also be written in the form 



^lu{x)^^''+^\x)dx- (\{x)dx = fA*i°VV 



x=l 
x=0 



(20) 



In particular, assume now that T-L is the n Hermite interpolant for tp in the unit interval, and 
apply (20) to both % and ip. Then the resulting right hand sides in the equations are equal, so 
that the left hand sides must also be equal. Thus, from (19), it follows that 

»i 







H{x)Ax^ / ip{x)Ax- / ^ik{x)^^^+^\x)dx 







(21) 







This can be scaled, to yield a formula that applies to the n*'* Hermite interpolant for 'ip on any 

1-D cell Xn "^^ X < Xn+l = Xn + h 



x„ + l 



Hix) Ax 



Xn + l 



ijj{x) dx — h 



fe+i 



h 



V'^'^+'^x) dx. 



(22) 



The following theorem is a direct consequence of this last formula 



Theorem 8. In one dimension, consider the computational domain fl = {x\a<x<b} and the 
grid: Xn — a + nh — where < n < N , h — ib — a)/N , and N > is an integer. Then, for any 
function tp G C'' , with integrahle (k + 1)*'* derivative 



{P,ltp){x)dx = / i;{x)dx 



i/''+'Hx)dx, 



where Ek = Ek{z) is the periodic function of period one defined by Ek{z) = /ife(zmodl). 

3.5.2. Stability for 1-D constant advection. Using Theorems [8] and [5] we can now prove stability 
for jet schemes in one dimension, and with a constant advection velocity. For simplicity we will 
also assume periodic boundary conditions. 

Theorem 9. Under the same hypothesis as in Theorem consider the constant coefficients 
advection equation (pt + v (px =0 in the computational domain Q,, with periodic boundary conditions 
(pib, t) = (p{a, t), and initial condition (f>{x, 0) = $(0). Approximate the solution by the jet scheme 
defined by the advect-and-project process of % 3.2 and the projection Pj^ , with a time step At 
satisfying At > 0(/i'^+^). Then this scheme is stable, in the sense described next. Let (pn be the 
numerical solution at time tn = n At, and consider a fixed time interval < i„ < T — for some 
fixed initial condition. Then: 



and 



lin bounded as h ^ 0, where < i < k. 



Notation: here \\-\\ is the norm in f2, and f^^^ denotes the i^'^ derivative of a function f. 
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Proof. For any Runge-Kutta scheme (and many other types of ODE solvers), the approximate 
advection operator is given by the shift operator 



{A(f>){x, t) = (j)[x ~vAt, t), 



(23) 



where we use the notation A = At+At. t — since At+^t, t is independent of time. Let the numerical 
solution at time t = tn = n At he denoted by Note that 0„ = (P^ o A)" $ G C'^. 

From Theorem [Hj Equation 23 and 0„+i = (P/l o A) 0„ it follows that 



< 



t+i) 



< 



(24) 



Furthermore, from Theorem [8] 

X((/.„+i) =X(0„) -/I'^+i 



El 



where Ai denotes the integral from a to h. From this, using the Cauchy-Schwarz inequality, we 
obtain the estimate 



\M{4>n+l)\<\M{4>n)\+h^+'C 



(25) 



where 



da: = (6- a) / {Ek{x)Y dx = 



b-a f{k + iy. 



n + 2 V(n+ 1)! 



From (25 1, using (24), it follows that 



|A^(0„)| < IMm + ^^h'^+'C 



$(fc+l) 



Now, for any of the (periodic) functions g 



< £ < k, we can write 



gix) = g+ I G{x-y)g^^\v)dy, 



(26) 



(27) 



where g is the average value for g and G is the Green's function defined by: G is periodic of period 
h- a, and G{x) = \ - ^ for < a; < 6 - a. Note that g = for 1 < £ < fc, and g = ^ M (</>„) 
for £ = 0. ° 



Now use (24), (27), and the Cauchy-Schwarz inequality, to find an /i-independent bound for 



— hence also for 



Repeat the process for (f)"n ^\ and so on — all the way down 



to 0!°'. In the last step, (|26|) is needed. 



□ 



Notice that Theorem [9] says nothing about the L°° norm of (p^n^^^ ■ A natural question is then: 
what can we say about (j)"n^^^ — beyond the statement in Equation (24 1? In particular, notice 
that there are stricter bounds on the growth of derivatives than the one given by Theorem [5], since 
!F is actually minimized in each cell individually. Can these, as well results similar to the one in 
Theorem [8) be exploited to obtain more detailed information about the behavior of the solutions 
(and derivatives) given by jet schemes? This is something that we plan to explore in future work. 



3.5.3. Stability for 1-D variable advection. An extension of this proof to the case of variable ad- 
vection will be considered in future work. 
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3.5.4. Stability for higher dimensions and constant coefficients. In several space dimensions, as- 
sume that the advection velocity does not involve rotation — specifically: the advected grid hyper- 
planes remain parallel to the coordinate hyperplanes. Then all the derivatives needed to compute 
the stability functional (as well as all the derivatives needed in the proof of Theorem [5]) remain 
defined after advection. Thus Theorem [5] can be used, as in the proof of Theorem[9) to control the 



growth of 
constant, 



where /? is as in Equation ([t]). In particular, when the advection velocity is 



< 



(28) 



follows. Just as in the one dimensional case, (28 1 alone is not enough to guarantee stability. For 



example: in the case with periodic boundary conditions, control over the growth of (f> still leaves 
the possibility of unchecked growth of partial means, i.e. components of the solution of the form 
4> = fj{x), where fj does not depend on the variable Xj. However, the Hermite interpolant for 
any such component reduces to the Hermite interpolant of a lower dimension, to which a lower 
dimensional version of the stability functional T applies — thus bounding their growth. The 
appropriate generalization of Theorem|8]to control the partial means of (p is the subject of current 
research. 



3.5.5. Stability for higher dimensions and variable coefficients. In the presence of rotation. The- 
orem [5] fails. Then, after advection, the derivatives involved are only defined in the sense of 



Definition 10 and the integrations by parts used to prove Theorem [5] are no longer justified. On 
the other hand (see the discussion at the beginning of § 2.4 1, the jumps in the derivatives that 
cause the failure are small. Hence we conjecture that, in this case, a "corrected" version of Theo- 
rem 5] will state that [H^] < J^[<j)] + C{h), where C{h) is a small correction that depends on the 
smallness of the jumps. In this case a simultaneous proof of stability and convergence might be 
possible — simultaneous because the "small jumps" property is valid only as long as the numerical 
solution is close to the actual solution. 



4. Numerical Results 

As a test for both the accuracy and the performance of jet schemes, we consider a version 
of the classical "vortex in a box" flow [31 [TT], adapted as follows. On the computational domain 
(x, y) € [0, 1]^, and for t e [0, tfinai], we consider the linear advection equation ([T]) with the velocity 
field 

(29) 



^, ( Trt\ ( sin^(7ra;) sin(27ru)\ 

u(a;, u, = cos (^) . ,\ ' . \, ^ {\ 
^ ' ^' ^ ^^-^ \-sin(27ra;) sin^(7ry)y 



This velocity field at t = is shown in Figure [TJ The maximum speed that ever occurs is 1. 
We prescribe periodic boundary conditions on all sides of the computational domain, and provide 



periodic initial conditions — note that the velocity field (29) is C°° everywhere when repeated 
periodically beyond [0, 1]^. This test is a mathematical analog of the famous "unmixing" exper- 
iment, presented by Heller [T^], and popularized by Taylor [57]. It models the passive advection 
of a solute concentration by an incompressible fiuid motion, on a time scale where diffusion can 
be neglected. The velocity field swirls the concentration forth and back (possibly multiple times) 
around the center point (i, i) in such a fashion that the equi-concentration contours of the solu- 
tion become highly elongated at maximum stretching. The parameters T and tfinai determine the 
amount of maximum deformation and the number of swirls. In all tests, we choose iflnai = iT, 
where £ d N. This implies that (t>{x, y, ifinai) = (t>{x, y, 0), i.e. at the end of each computation, the 
solution returns to its initial state. Below, we compare the accuracy and convergence rates of jet 



schemes with WENO methods (§ 4.1), investigate the accuracy and stability of different versions 



of jet schemes (§ 4.2), demonstrate the performance of jet schemes on a benchmark test (§ 4.3), 



and finally compare the computational cost and efficiency of all the numerical schemes considered 



(§ 4.4) 
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0.2 0.4 0.6 0.8 1 



Figure 1. Streamlines and quiver plot of the velocity field used for the test cases. 



4.1. Test of Accuracy and Convergence Rate. We first test the relative accuracy and conver- 
gence rate of jet schemes in comparison with classical WENO approaches. For this test, we choose 
ifinai = T = I and initial conditions (f>{x, y, 0) = cos(27rx) cos(47r?/). This choice of parameters 
yields a moderate deformation, and the smallest structures in the solution are well resolved for 
h < Vso- In all error convergence tests, the error with respect to the true solution at tfinai is 
measured in the L°° norm. 

The results of the numerical error analysis are shown in Figure 



2l The jet schemes considered 
here are the ones based on the projection , given in Definition 11 Specifically, we consider a 
bi-linear scheme (fc = 0), denoted by black dots; a bi-cubic schenie(fc = 1), implemented using 
analytical differentiation for the partial 1-jet, as described in § |3.3.1[ marked by black squares; 
and the bi-quintic scheme {k = 2) described in the second bullet point in Example [5] denoted by 
black triangles. 

As reference schemes, we consider the classical WENO finite difference schemes described in 
|15) . We use schemes of orders 1, 3, and 5, constructed as follows. Unless otherwise noted, time 
stepping is done with the maximum time step that stability admits. The first order version is the 
simple 2-D upwind scheme, using forward Euler in time (with a time step At = ^^h). The time 
stepping of the third order WENO is done using the Shu-Osher scheme |1S] (with At = h). Lacking 
a simple fifth order SSP Runge-Kutta scheme [H [TU] , the fifth order WENO is advanced forward 
in time using the same third order Shu-Osher scheme, however, with a time step At — hi . This 
yields a globally order 0{h^) scheme. For both WEN03 and WEN05, the parameters e = 10~^ 
and p — 2 {in the notation of |15j) are used. These choices are commonly employed when WENO 
is used in a black-box fashion (i.e. without adapting the parameters to the specific solution). In 
Figure [2j the numerical errors incurred with the WENO schemes are shown by gray curves and 
symbols, namely: dots for upwind; squares for WEN03; and triangles for WEN05. 

This test demonstrates the potential of jet schemes in a remarkable fashion. While the first 
order jet scheme is very similar to the 2-D upwind method, the jet schemes of orders 3 and 5 are 
strikingly more accurate than their WENO counterparts. For equal resolution h, the errors for 
the bi-cubic jet scheme are about 100 times smaller than with WEN03, and for the bi-quintic 
jet scheme, the errors are about 1000 times smaller than with WEN05. This implies that a high 
order jet scheme achieves the same accuracy as a WENO scheme (of the same order) with a 
resolution that is about 4 times as coarse. This reflects the fact that jet schemes possess subgrid 
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resolution, i.e. structures of size less than h can be represented due to the high degree polynomial 
interpolation pT| . 

4.2. Accuracy and Stability of Grid-Based Finite Difference Schemes. In a second test, 
the accuracy and stability of schemes that use grid-based finite difference reconstructions of the 



higher derivatives are tested. As in the test described in § 4.1 we choose the deformation T — 1, 
however, now ifinai — 20, i.e. the solution is moved back and forth multiple times. Through this 
approach, enough time steps are taken so that unstable approaches in fact show their instabilities. 

We compare jet schemes of order 3 (fc = 1) and order 5 (fc = 2), each in two versions: first, 
we consider the approaches based on the projection (see Definition 111, to which the stability 



arguments in § |3.5| apply. Second, we implement schemes that construct the partial /c-jet from 
t he to tal fc-jet, using optimally local grid-based finite difference approximations. As described in 
these types of schemes are generally less costly than approaches based on P^^ , and therefore 



2.5 



wort 



1 investigating (see also § 4.4) 



Specifically, we consider a bi-cubic scheme that tracks the total 1-jet, and uses the grid-based 
finite difference approximation of (f)xy described in Example [3j Furthermore, we implement a 
bi-quintic scheme that tracks the total 2-jet, and approximates (l)xxyi i'xyyj and 4>xxyy using grid- 
based finite differences, as follows. In a cell Q — [0, /i]^, with the index q S {0, 1}^ corresponding 
to the vertex ^ qh, the derivatives at the vertex (0, 0) are approximated by: 



h \ ^xy ^ ^xy J ' 2h y ^vv ' ^yy ^yy ' ^yy 

It can be verified by Taylor expansion that these are 0{h^^'') accurate approximations to the 
respective derivatives of order s. The corresponding approximations at the vertices (1,0), (0, 1) 
and (1,1) are obtained by symmetry. 

The convergence errors for these jet schemes are shown in Figure [3] As already observed in 



§ 4.1 the bi-cubic (shown by black squares) and bi-quintic (denoted by black triangles) jet schemes 
that are based on tracking the partial A:-jet are stable and {2k + 1)*** order accurate. The bi-cubic 
jet scheme that approximates (pxy by grid-based finite differences (indicated by open squares) 
is stable and third order accurate as well, albeit with a larger error constant than the "pure" 
version based on . In contrast, the bi-quintic jet scheme that approximates (pxxy, 4'xyy: and 
(t^xxyy by grid-based finite differences (denoted by open triangles) is unstable. Interestingly, the 
instability is extremely mild: even though many time steps are taken in this test, the instability 
only shows up for h < 0.008. This phenomenon is both alarming and promising. Alarming, 
because it demonstrates that with jet schemes that are not based on the partial fc-jet, stability 
is not assured, and generalizations of the stability results given in § |3.5| are needed. Promising, 
because it is plausible that these jet schemes could be stabilized without significantly having to 
diminish their accuracy. 

4.3. Test of Performance on a Level-Set-Type Example. In this test we assess the practical 
accuracy of the considered numerical approaches, by following a curve of equi-concentration of the 
solution. The parameters are now ifinai = T = 10, and the initial conditions are given by a periodic 
Gaussian hump: (j){x, y, 0) =J2i,jez9ix-i, y-j), where g{x, y) = exp(-10((a;-Xo)^-K(2/-yo)^)) 
with a;o — 0.5 and yo = 0.75. We examine the time-evolution of the contour T{t) — {{x, y) G 
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Figure 2. Numerical convergence rates for jet schemes of 
orders 1, 3, and 5, in comparison with WENO schemes of 
the same orders. 
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Figure 3. Numerical convergence rates for jet schemes of 
orders 3 and 5, comparing schemes based on the full partial 
jet with schemes using grid-based finite differences. The 
latter type of approach turns out unstable for the order 5 
jet scheme. 
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Figure 4. Swirl test: velocity field and 
initial conditions. 
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Figure 5. Swirl test: (f) = (f)c 
maximum deformation {t = 



contour at 
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Figure 6. Swirl test: 
the final time (t^T). 



= 0c contour at 
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ri : (j){x, y, t) = (j)c} that corresponds to the concentration (j>c = exp(— lOr^), where r = 0.15. 
The initial (and final) contour r(0) is almost an exact circle of radius r, centered at (xq, yo). At 
maximum deformation, r(-|^) is highly elongated. 

This test is well-known in the area of level set approaches [22]. For these, only one specific 
contour is of interest, and it is common practice to modify the other contours, e.g. by adding a 
reinitialization equation |26| . or by modifying the velocity field away from the contour of interest, 
using extension velocities T. Since here the interest lies on more general advection problems 
([T]), such as the convection of concentration fields, no level-set-method-specific modifications are 
considered. However, this does not mean that these procedures cannot be applied in the context 
of jet schemes. In fact, proper combinations of jet schemes with reinitialization are the subject of 
current research. 

We apply the bi-cubic and the bi-quintic jet schemes, as well as the WEN05 scheme (all as 



described in § 4.1|, on a grid of resolution h — ^/ loo- The resulting equi-concentration contours are 
shown in Figure |4| (initial conditions), Figurejs] (maximum deformation), and Figure [6] (final state). 
The contour obtained with WEN05 is thin and black, the contour obtained with the bi-cubic jet 
scheme is medium thick and blue, and the contour obtained with the bi-quintic jet scheme is thick 
and red. The true solution (approximated with high accuracy using Lagrangian markers), is shown 
as a gray patch. For the considered resolution, both jet schemes yield more accurate results than 
WEN05. In fact, the fifth order jet scheme yields an almost flawless approximation to the true 
solution, on the scale of interest. As alluded to in § |4.2[ the subgrid resolution of the jet schemes 



is of great benefit in resolving the thin elongated structure. 

4.4. Computational Cost and EfRciency. In order to get an impression of the relative com- 
putational cost and efficiency of the presented schemes, we measure the CPU times that various 
versions of jet schemes and WENO require to perform given tasks. We use the same test as in 
and apply the following versions of jet schemes: (a) jet schemes of orders 1, 3, and 5, that 



4.1 



are based on the partial fc-jet; and (b) jet schemes of orders 3 and 5, that are based on the total 
fc-jet and use grid-based finite difference approximations for the missing derivatives of the partial 



fc-jet (see § 4.2 for a description of these schemes). In addition, we consider the following versions 
of WENO and finite difference schemes to serve as reference methods: simple linear 2-D upwind- 
ing with forward Euler in time (with At = ':^^)^ WEN03 [15] advanced with the third order 
Shu-Osher method [IS], and WEN05 with the Cash-Karp method [S] used for the time stepping. 



Note that unlike the case where convergence is tested (see § 4.1), for a fair comparison of CPU 
times a true fifth order time stepping scheme must be used. While the Cash-Karp method is not 
SSP, for the considered test case no visible oscillations are observed. Furthermore, in order to give 
a fair treatment to WENO, one must be careful with the choice of the parameter e that WENO 
schemes require [H]. The meaning of the parameter is that for e ^ /i^, spurious oscillations are 
minimized at the expense of a degradation of convergence rate, while for e ^ h? , the schemes 
reduce to (unlimited) linear finite differences of the respective order. We therefore consider two 
extremal cases for both WEN03 and WEN05: one with a tiny e — 10"^", and one that is linear 
finite differences (called "linear FD3" and "linear FD5") with upwinding (with all limiting routines 
removed). Practical WENO implementations will be somewhere in between these two extremal 
cases. 

As a first comparison of the computational costs of these numerical approaches, we run all of 
them for a fixed resolution [h = Vi5o), and record the resulting L°° errors and CPU times, as 
shown in Table [l] One can see a clear separation between methods of different orders both in 
terms of accuracy and in terms of CPU times. When comparing approaches of the same order one 
can observe that for identical resolutions, jet schemes incur a significantly larger computational 
effort than WENO and linear finite difference schemes. However, this increased cost comes at 
the benefit of a significant increase in accuracy. Similarly, it is apparent that the use of grid- 
based finite difference approximations for parts of the fc-jet can reduce the cost of jet schemes 
significantly. However, this tends to come at an expense in accuracy, in addition to the fact that 
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Order 


Approach 


L°° error 


CPU time/sec. 


1 


bi-linear jet scheme 


1.69 X 10-^ 


0.400 


1 


upwind 


1.92 X 10-1 


0.456 


3 


bi-cubic jet scheme 


1.35 X 10"'' 


7.91 


3 


bi-cubic jet scheme grid-based FD 


2.31 X 10"'' 


5.99 


3 


WEN03/RK3 (e = 10-1°) 


1.21 X 10^2 


1.61 


3 


Unear FD3/RK3 


1.54 X 10^3 


1.15 


5 


bi-quintic jet scheme 


8.23 X 10-** 


103 


5 


bi-quintic jet scheme grid-based FD 


1.32 X 10-^ 


36.3 


5 


WEN05/RK5 (e = IQ-^") 


1.25 X 10--* 


19.3 


5 


Hnear FD5/RK5 


2.15 X 10-^ 


19.1 



Table 1. CPU times for a fixed resolution h — Yiso with various jet schemes and WENO schemes. 



Order 


Approach 


Res. h 


L°° error 


time/sec. 


1 


bi-hnear jet scheme 


V 1600 


1.89 X 10-2 


441 (*) 


1 


upwind 


V 1600 


2.28 X 10-2 (*) 


482 (*) 


3 


bi-cubic jet scheme 


V 167 


9.86 X 10"^ 


11.4 


3 


bi-cubic jet scheme grid-based FD 


V 200 


9.86 X 10-5 


13.0 


3 


WEN03/RK3 (e = lO-^°) 


V 1350 


9.61 X 10-5 


1110 


3 


linear FD3/RK3 


V38O 


9.65 X 10-5 


18.7 


5 


bi-quintic jet scheme 


V 35 


9.76 X 10-5 


1.55 


5 


bi-quintic jet scheme grid-based FD 


V40 


8.62 X 10-5 


0.821 


5 


WEN05/RK5 (e = 10-i°) 


V 158 


9.71 X 10-5 


21.4 


5 


Hnear FD5/RK5 


V 110 


9.87 X 10-5 


11.7 



Table 2. CPU times required to achieve an L°° error of less than 10 *. Note the lower accuracy of the 
first order schemes. 



the resulting scheme could become unstable, as demonstrated in § |4.2[ In fact, the CPU times for 
the bi-quintic jet scheme with grid-based finite differences shown in Table [T] must be interpreted 
as a guideline only, since without further stabilization this version of the method is of no practical 
use. 

As a second test we compare the schemes in terms of their actual efficiency. To that end, we 
measure the CPU times that the schemes require to achieve a given accuracy, i.e. for each scheme 
we choose the coarsest resolution such that the L°° error is below 10"* , and for this resolution 
we report the required CPU time. The results of this test are shown in Table [2] Note that the 
first order schemes were computed at a maximum resolution oih = Vieoo- A simple extrapolation 
reveals that in order to reach an accuracy of 10-"*, CPU times of more than 50 years would be 
required. In contrast, the third and fifth order methods yield reasonable CPU times, except for 
WEN03 with e = 10-^°, which due to the degraded convergence rate requires an unreasonably 
fine resolution to reach the target accuracy. In fact, for this very smooth test case, linear finite 
differences present themselves as rather efficient methods, which achieve descent accuracies at 
a very low cost per time step. However, the results show that jet schemes are actually more 
efficient than their WENO/finite difference counterparts. Their higher cost per time step is more 
than outweighed by the significant increase in accuracy. Specifically, third order jet schemes are 
observed to be more efficient than finite differences by a factor of 1.5, and fifth order jet schemes 
by a factor of at least 7.5. 

The computational efficiency of jet schemes goes in addition to their structural advantages, such 
as optimal locality and subgrid resolution. A more detailed efficiency comparison of jet schemes 
with WENO, and also with Discontinuous Galerkin approaches, is given in 
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5. Conclusions and Outlook 

The jet schemes introduced in this paper form a new class of numerical methods for the linear 
advection equation. They arise from a conceptually straightforward formalism of advect-and- 
project in function spaces, and provide a systematic methodology for the construction of high 
(arbitrary) order numerical schemes. They arc based on two ingredients: an ODE solver to track 
characteristic curves, and Hermite interpolation. Unlike finite volume methods or discontinuous 
Galerkin approaches, jet schemes have (currently) a limited area of application: linear advection 
equations that are approximated on rectangular grids. However, under these circumstances, they 
have several interesting advantages: a natural treatment of the flow nature of the equation and 
of boundary conditions, no CFL restrictions, no requirement for SSP ODE solvers, and optimal 
locality, i.e. the update rule for the data on a grid point uses information in only a single grid cell. 

Jet schemes achieve high order by carrying a portion of the jet of the solution as data. To 
derive an update rule in time for the data, the notion of superconsistency is introduced. This 
concept provides a way to systematically inherit update procedures for the derivatives from an 
approximation scheme (ODE solver) for the characteristic curves. Various alternative approaches 
are presented for how to track different portions of the jet of the solution, and for how to specifi- 
cally implement superconsistent schemes. One example is the grid-based reconstruction of higher 
derivatives, which on the one hand can dramatically reduce the computational cost, but on the 
other hand could destabilize the scheme. Another example is the idea of e-finite differences, which 
requires diligent considerations of round-off errors, but can lead to simpler implementations and 
also reduce the computational cost. 

The formal equivalence of superconsistent schemes to advect-and-project approaches in func- 
tion spaces gives rise to a notion of stability: the projection step (which is based on a piece-wise 
application of Hermite interpolation) minimizes a stability functional, while the increase by the 
advection step can be appropriately bounded. The stability functional gives control over certain 
derivatives of the numerical solution, and thus bounds the occurrence of oscillations. In this paper, 
a full proof of stability for the 1-D constant coefficient case is provided. 

Investigations of the numerical error convergence and performance of high order jet schemes 
show that they tend to be more costly, but also strikingly more accurate than classical WENO 
schemes of the same respective orders. It is observed that, for the same resolution, a jet scheme 
incurs 2 7 times the computational cost of a WENO scheme. However, the numerical experiments 
also show that a jet scheme achieves the same accuracy as a WENO scheme (of the same order) at 
a 3-4 times coarser grid resolution. Efficiency comparisons reveal that jet schemes tend to achieve 
the same accuracy as an efficient WENO scheme at a 1.5 7 times smaller computational cost. 
In comparison with WENO schemes with a sub-optimal parameter choice, the efficiency gains of 
jet schemes are factors of 14-85. An application to a model of the passive advection of contours 
demonstrates the practical accuracy of jet schemes. A particular fcatiue of jet schemes is their 
potential to represent structures that are smaller than the grid resolution. 

This paper should be seen as a first step into the class of jet schemes. Their general interpreta- 
tion as advect and project approaches, their conceptual simplicity, and their apparent accuracy 
make them promising candidates for numerical methods for certain important types of advection 
problems. However, many questions remain to be answered, and many aspects remain to be in- 
vestigated, before one could consider them generally "competitive". The listing below gives a few 
examples of some of the aspects where further research and development is still needed. 

Analysis. A general proof of stability in arbitrary space dimensions, and for variable coeffi- 
cients, is one important goal for jet schemes. A fundamental understanding of the stability of these 
methods is particularly important when combined with grid-based finite difference approximations 
to higher derivatives. 
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Efficiency. Some ideas presented liere (e.g. grid-based or e-finite differences) lead to significant 
reductions in implementation effort and computational cost, but they also indicate that the sta- 
bility of such new approaches is not automatically guaranteed. One can safely expect that there 
arc many other ways to implement jet schemes more easily and more efficiently. The investigation 
and analysis of such ideas is an important target for future research. 

Applicability. The field of linear advection problems on rectangular grids is not as limited as 

it may seem. Many examples of interface tracking fall into this area. Of particular interest is the 
combination of jet schemes with adaptive mesh refinement techniques. The optimal locality of jet 
schemes promises a straightforward generalization to adaptive meshes. 

Generality. Since all that jet schemes need are a method for tracking the characteristics, and 

an interpolation procedure, they are not fundamentally limited to rectangular grids. Similarly, 
they could be generalized to allow for non-conforming boundaries, and thus apply to general 
unstructured domains. Furthermore, it is plausible that the idea of tracking derivatives can be 
advantageous for nonlinear Hamilton- Jacobi equations, hyperbolic conservation laws, or diffusion 
problems as well. In all these cases new challenges arise, such as the occurrence of shocks or the 
absence of characteristics. 
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